Worflow

runZinb <- T
runClus <- T
NCORES <- 10

EAP: if we are going to ultimately use these functions, it would be good to add them to the package for zinbwave.

FP: yes, I agree.

compute.zinb.loglik <- function(Y, zinb){
  mu = t(getMu(zinb)) 
  theta = getTheta(zinb)
  theta_mat = matrix(rep(theta, ncol(Y), ncol = ncol(Y)))
  pi = t(getPi(zinb))
  log( pi * (Y == 0) + (1 - pi) * dnbinom(Y, size = theta, mu = mu) )
}

compute.deviance.residuals <- function(Y, zinb){
  mu_hat = t(getMu(zinb)) 
  pi_hat = t(getPi(zinb)) 
  Y_hat = (1 - pi_hat) * mu_hat
  ll = compute.zinb.loglik(Y, zinb)
  sign = 1*(Y - Y_hat > 0)
  sign[sign == 0] = -1
  sign * sqrt(-2 * ll)
}

TL;DR

This document is a temptative of workflow with the following steps

• Dimensionality reduction using zinbwave to get W which should capture the biology,
• Cluster cells using clusterExperiment on W to get the cluster labels,
  
• Get lineage using slingshot on W and cluster labels from clusterExperiment,
  
• Get DE genes between lineages/clusters.

Along the worflow, use deviance residuals as adjusted values.

EAP: small request. Can everyone put a line between the beginning of a r chunk and text? It makes it nicely formated for my text editor.

Steps of the workflow

Create a SummarizedExperiment object

Along the workflow, we want to use a unique SummarizedExperiement object carrying all the data we need.

EAP: I have updated the code to pull from a dataset on the repos that is created with the createData.R file. For now, I am filtering to the top 1000 most variable genes there, though we might want to add that to the code for the article. This will be slightly different data from Russell’s, which didn’t use all of the samples. We can adjust that decision later, or just compare the samples that are the same. (Russell’s clusterLabels are in the meta data)

EAP: zinbFit doesn’t accept data.frame objects, so currently have to have a data.matrix command. Should it be changed so that it does?

#counts<-read.table("../data/oeCufflinkCountData.txt",sep="\t",header=TRUE)
core <- read.table("../data/oeCufflinkCountData_1000Var.txt",
                   sep = "\t", header = TRUE)
core <- data.matrix(core)
metadata <- read.table("../data/oeMetadata.txt", sep = "\t", header = TRUE)
# symbol for samples missing from original clustering
metadata$clusterLabels[is.na(metadata$clusterLabels)] <- -2

Here we only look at the 1000 most variable genes. EAP: see note above, I’ve commented out the filtering and added it to the createData.R for now.

batch <- metadata$Batch

Cells have been processed in 18 different batches

col_batch = rep(brewer.pal(9, "Set1"), 2)
names(col_batch) = unique(batch)
table(batch)

## batch
## GBC08A GBC08B GBC09A GBC09B    P01    P02   P03A   P03B    P04    P05 
##     41     47     43     38     34     49     73     52     24     23 
##    P06    P10    P11    P12    P13    P14    Y01    Y04 
##     53     51     54     51     60     49     65     42

We have qc measures from the data

qc <- metadata[, !names(metadata) %in% c("Batch", "Experiment", "clusterLabels")]
head(qc, 2)

##                  NREADS NALIGNED  RALIGN TOTAL_DUP    PRIMER
## OEP01_N706_S501 3313260  3167600 95.6035   47.9943 0.0154566
## OEP01_N701_S501 2902430  2757790 95.0167   45.0150 0.0182066
##                 PCT_RIBOSOMAL_BASES PCT_CODING_BASES PCT_UTR_BASES
## OEP01_N706_S501               2e-06         0.200130      0.230654
## OEP01_N701_S501               0e+00         0.182461      0.201810
##                 PCT_INTRONIC_BASES PCT_INTERGENIC_BASES PCT_MRNA_BASES
## OEP01_N706_S501           0.404205             0.165009       0.430784
## OEP01_N701_S501           0.465702             0.150027       0.384271
##                 MEDIAN_CV_COVERAGE MEDIAN_5PRIME_BIAS MEDIAN_3PRIME_BIAS
## OEP01_N706_S501           0.843857           0.061028           0.521079
## OEP01_N701_S501           0.914370           0.033350           0.373993
##                 CreER ERCC_reads
## OEP01_N706_S501     1      10516
## OEP01_N701_S501  3022       9331

clus.labels <- metadata[, "clusterLabels"]

In published work (FP: add ref), cells have been clustered into 14 different clusters

col_clus <- c("transparent", brewer.pal(12, "Set3"), brewer.pal(8, "Set2"))
col_clus <- col_clus[1:length(unique(clus.labels))]
names(col_clus) <- sort(unique(clus.labels))
table(clus.labels)

## clus.labels
##  -2   1   2   3   4   5   7   8   9  10  11  12  14  15 
## 233  91  25  56  40  96  60  28  79  26  22  35  26  32

Batches are kind of confounded with the biology

table(data.frame(batch = as.vector(batch),
                 cluster = clus.labels))

##         cluster
## batch    -2  1  2  3  4  5  7  8  9 10 11 12 14 15
##   GBC08A  5  0  2 12  9  0  0  0  0  0  2  0  2  9
##   GBC08B 14  0  7  5  3  0  0  0  1  2  4  0  5  6
##   GBC09A 13  0  1  5  9  0  0  0  1  1  0  0  6  7
##   GBC09B 21  0  2  2  7  0  0  0  3  0  0  0  3  0
##   P01     9  0  2  4  3 15  1  0  0  0  0  0  0  0
##   P02     6  2  0  9  3 15  3  3  2  3  0  2  1  0
##   P03A   36  3  0  3  0 12  2  9  4  2  0  2  0  0
##   P03B   19  1  2  1  1 11  1  2 10  1  1  2  0  0
##   P04    10  0  0  0  0 11  1  0  1  1  0  0  0  0
##   P05     3  0  0  0  1 11  3  0  1  0  2  2  0  0
##   P06    14  1  2  3  0  8  2  4  8  4  1  2  2  2
##   P10    15  3  1  4  0  4  5  9  2  0  2  5  0  1
##   P11    10  2  1  1  0  1  5  1 22  3  1  6  0  1
##   P12    11  0  2  0  0  4 10  0  8  2  3  6  4  1
##   P13    13  1  2  4  0  4 15  0  4  5  6  1  3  2
##   P14    10  0  0  1  2  0 12  0 12  2  0  7  0  3
##   Y01    14 47  1  1  2  0  0  0  0  0  0  0  0  0
##   Y04    10 31  0  1  0  0  0  0  0  0  0  0  0  0

We have 849 cells.

dim(core)

## [1] 1000  849

core[1:3, 1:3]

##        OEP01_N706_S501 OEP01_N701_S501 OEP01_N707_S507
## Cbr2              5799            3638            1448
## Cyp2f2            2158            2027            1078
## Gstm1             8763            7221            3581

Let’s create a SummarizedExperiment object to store the raw counts and information about the data, that is batches, published labels, and quality control measures.

se <- SummarizedExperiment(assays = list(rawCounts = core),
                           colData = metadata)

Dimensionality reduction adjusting for gene and cell level covariates

To cluster and get lineages we want to reduce the dimension of the data. We are going to use zinbwave to do so. First, let’s fit zinbwave with first K = 0 to compute normalized values (i.e. deviance residuals) adjusted for batches. We could also adjust for gene length or GC content here. We then fit zinbwave to get the dimensionality reduced matrix W specifying the number of dimension K = 50. Eventually, we will call zinbwave just once where we would have an argument in zinbFit like “compute_normalized_values” in c(TRUE, FALSE). For K = 0 and K = 50, we correct for batch effect including batches in X.

fn0 <- '../data/zinb_k0_batch.rda'
if (runZinb & !file.exists(fn0)){
  mod <- model.matrix( ~ batch)
  print(system.time(zinb0 <- zinbFit(core, ncores = NCORES,
                                     K = 0, X = mod)))
  save(zinb0, file = fn0)
}else{
  load(fn0)
}

fn50 <- '../data/zinb_k50_batch.rda'
if (runZinb & !file.exists(fn50)){
  mod <- model.matrix( ~ batch)
  print(system.time(zinb50 <- zinbFit(core, ncores = NCORES,
                                      K = 50, X = mod)))
  save(zinb50, file = fn50)
}else{
  load(fn50)
}

Normalized values

We use deviance residuals as normalized values for visualization. FP: explain rational: K=0 so residuals capture the bio adjusting for batch. Let’s check that deviance residuals look ok.

res <- compute.deviance.residuals(core, zinb0)
res[1:3,1:3]

##        OEP01_N706_S501 OEP01_N701_S501 OEP01_N707_S507
## Cbr2          4.541942       -4.431585       -4.239279
## Cyp2f2        4.316532        4.301158       -4.146744
## Gstm1         4.629884        4.583829       -4.429794

Boxplot of the normalized values for each cell. It seems that correction for batches is ok.

res_order <- res[, order(as.numeric(batch))]
col_order <- as.numeric(batch)[order(as.numeric(batch))]
boxplot(res_order, main='Boxplot of normalized values\ncolor=batch',
        col = col_order, staplewex = 0, outline = 0, border = col_order, xaxt = 'n')

PCA on the normalized values where color are for batches on the left and previously found clusters on the right. We want no clustering on the left side and clustering on the right side.

pca <- prcomp(t(res))
par(mfrow = c(1,2))
plot(pca$x, col = col_batch[batch], pch = 20,
     main="PCA of normalized values\ncolor=batch")
plot(pca$x, col = col_clus[as.character(clus.labels)], pch = 20,
     main = "PCA of normalized values\ncolor=cluster")

par(mfrow = c(1,1))

Let’s add the normalized values as an assay dataset in our SummarizedExperiment object.

assays(se)[['normalizedValues']] <- res

W

Let’s check that performing MDS on W we have something coherent with published clusters.

W <- getW(zinb50)
d <- dist(W)
fit <- cmdscale(d, eig = TRUE, k = 2)
plot(fit$points, col = col_clus[as.character(clus.labels)], main = 'MDS', pch = 20,
     xlab = 'Component 1', ylab = 'Component 2')
legend(x = 'bottomright', legend = unique(names(col_clus)), cex = .5,
       fill = unique(col_clus), title = 'Sample')

Let’s add W to the colData of our SummarizedExperiment object.

W <- data.frame(W)
colnames(W) <- paste0('W', 1:ncol(W))
colData(se) <- cbind(colData(se), W)

FP: what do you think of implemented a “compute_normalized_values” argument to zinbFit for summarizedExperiment objects where zinbFit would take as input a SummarizedExperiment object and return a summarizedExperiment object with slots added for normalized values and W?

Cluster of cells

We use clusterExperiment with W.

EAP: Fanny, we should really be using the RSEC function for all the steps of the clustering. I’ve put in the code, but I set eval=FALSE because I don’t have the data to test it. I also am not sure that we want to set the options for subsampleArgs and seqArgs since those are suppose to be sort of esoteric concerns…

FP: What about the following call to RSEC?

W <- colData(se)[, grepl('^W', colnames(colData(se)))]
W <- as.matrix(W)
fn <- '../data/RSEC_W.rda'
if (runClus & !file.exists(fn)){
  print(system.time(ceObj <- RSEC(t(W), k0s = 4:15, 
                                  alphas = c(0.3), betas = c(0.9),
                                  clusterFunction = "hierarchical01",
                                  minSizes = 5,
                                  ncores = NCORES,
                                  combineProportion = 0.7,
                                  mergeMethod = "locfdr",
                                  mergeCutoff = 0.01)))
  save(ceObj, file = fn)
}else{
  load(fn)
}

## Note: Merging will be done on ' combineMany ', with clustering index 1

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

##     user   system  elapsed 
## 4823.683    7.293 1038.270

plotClusters(ceObj, colPalette = c(bigPalette, rainbow(100)))

plotCoClustering(ceObj)

## Warning in .makeColors(clusters, colors = bigPalette): too many clusters to
## have unique color assignments

table(primaryClusterNamed(ceObj))

## 
##  -1  m1 m10 m11 m12 m13 m14 m15 m16  m2  m3  m4  m5  m6  m7  m8  m9 
## 520  56  10  58  12   5   6   7   5  13   9   5  70  39  16  13   5

sum(primaryCluster(ceObj) == -1)

## [1] 520

FP: Elizabeth, we are working with the W here, does the locfdr make sense in this context?

#re-does merging simpling to make plot 
mergeClusters(ceObj, mergeMethod = "locfdr",
              plotInfo = "mergeMethod", cutoff = 0.05)

## Note: Merging will be done on ' combineMany ', with clustering index 2

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning: glm.fit: fitted rates numerically 0 occurred

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

## Warning in locfdr::locfdr(tstats, plot = 0): CM estimation failed, middle
## of histogram non-normal

plotHeatmap(ceObj, clusterSamplesData = "dendrogramValue", 
            breaks = .99)

plot(fit$points, col = col_clus[as.character(clus.labels)],
     main = 'MDS W, color = published clusters', pch = 20,
     xlab = 'Component1', ylab = 'Component2')

palDF <- ceObj@clusterLegend[[1]]
pal <- palDF[, 'color']
names(pal) <- palDF[, 'name']
plot(fit$points, col = pal[primaryClusterNamed(ceObj)],
     main = 'MDS W, color = our clusters', pch = 20,
     xlab = 'Component1', ylab = 'Component2')

table(data.frame(published = clus.labels,
                 ours = primaryCluster(ceObj)))

##          ours
## published -1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
##        -2 83 24  1 16 51  2 22  3  1  0  2 12  5  3  0  0  3  2  0  1  2
##        1  25 54  1  3  3  4  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##        2   2  0  0  0  0  0  0  0  0 20  0  0  0  0  0  0  0  3  0  0  0
##        3  10  0  0  4  0  0  0  0 37  1  0  0  0  0  0  3  0  0  0  0  1
##        4   9  0  0  0 28  0  2  0  0  0  1  0  0  0  0  0  0  0  0  0  0
##        5  28 18  3  2 38  1  0  2  0  0  2  2  0  0  0  0  0  0  0  0  0
##        7  10  0  0  0 47  0  1  0  0  0  0  1  0  0  0  0  1  0  0  0  0
##        8  14 12  0  1  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##        9  40  2  0  0  3  0  1  0  0  0  0  2 28  0  3  0  0  0  0  0  0
##        10 19  0  0  0  0  0  0  0  0  0  0  0  4  0  2  0  1  0  0  0  0
##        11  9  0  0  0  0  0  0  0  6  0  0  0  0  2  0  3  0  0  0  0  2
##        12  4  0  0  0  0  0  0  0  0  0  0  0 31  0  0  0  0  0  0  0  0
##        14 16  1  0  0  1  0  0  0  0  0  0  0  0  1  1  0  0  2  0  4  0
##        15 18  0  0  2  0  0  0  0  0  0  0  0  0  0  0  0  0  0 12  0  0

So, let’s look at a heatmap on normalized values.

sampleData <- data.frame(ours = primaryCluster(ceObj))

plotHeatmap(assays(se)$normalizedValues,
            main = 'Normalized values, 1000 most variable genes',
            clusterSamplesData = ceObj@dendro_samples,
            sampleData = as.matrix(sampleData))

## Warning in .Method(..., deparse.level = deparse.level): number of rows of
## result is not a multiple of vector length (arg 1)

Get lineages

The goal of this section is to see if we need to refit zinbwave when we want to run slingshot. We run slingshot in the supervised and unsupervised mode and we try k=2, k=3, k=4 in W.

FP: We need to work on the cluster coloring (!) and the matching with the clusters you found for the paper with Russell, maybe by changing the parameters in clusterMany?

Kvec <- c(2, 3, 4)

Use previous W, no re-fitting

Unsupervised

our_cl <- as.factor(primaryClusterNamed(ceObj))
cl = our_cl[our_cl != "-1"]
for (k in Kvec){
  X <- W[our_cl != "-1", 1:k]

  lineages <- get_lineages(X, clus.labels = cl, start.clus = "m1")
  curves <- get_curves(X, clus.labels = cl, lineages = lineages)
  plot_curves(X, cl, curves, col.clus = pal)
  plot_tree(X, cl, lineages, col.clus = pal)

  print(paste0("K=", k))
  print(lineages$lineage1)
  print(lineages$lineage2)
  print(lineages$lineage3)
  print(lineages$lineage4)
}

## [1] "K=2"
## [1] "m1"  "m5"  "m16" "m11" "m13" "m17" "m9" 
## [1] "m1"  "m5"  "m16" "m11" "m13" "m17" "m19"
## [1] "m1"  "m5"  "m4"  "m2"  "m10" "m7" 
## [1] "m1"  "m3"  "m20" "m15" "m8"

## [1] "K=3"
##  [1] "m1"  "m5"  "m2"  "m7"  "m10" "m4"  "m11" "m16" "m3"  "m20" "m15"
## [12] "m8" 
##  [1] "m1"  "m5"  "m2"  "m7"  "m10" "m4"  "m11" "m13" "m17" "m9" 
##  [1] "m1"  "m5"  "m2"  "m7"  "m10" "m4"  "m11" "m16" "m12" "m14"
##  [1] "m1"  "m5"  "m2"  "m7"  "m10" "m4"  "m11" "m13" "m17" "m19"

## [1] "K=4"
##  [1] "m1"  "m5"  "m2"  "m10" "m4"  "m16" "m11" "m13" "m17" "m9" 
##  [1] "m1"  "m5"  "m2"  "m10" "m4"  "m16" "m11" "m13" "m17" "m19"
## [1] "m1"  "m5"  "m2"  "m10" "m4"  "m16" "m11" "m6" 
## [1] "m1"  "m5"  "m2"  "m10" "m4"  "m16" "m12" "m14"

Supervised

for (k in Kvec){
  X <- W[our_cl != "-1", 1:k]

  lineages <- get_lineages(X, clus.labels = cl, start.clus = "m1",
                           end.clus = c("m8", "m4"))
  curves <- get_curves(X, clus.labels = cl, lineages = lineages)
  plot_curves(X, cl, curves, col.clus = pal)
  plot_tree(X, cl, lineages, col.clus = pal)

  print(paste0("K=", k))
  print(lineages$lineage1)
  print(lineages$lineage2)
  print(lineages$lineage3)
  print(lineages$lineage4)
}

## [1] "K=2"
## [1] "m1"  "m5"  "m16" "m11" "m13" "m17" "m9" 
## [1] "m1"  "m5"  "m16" "m11" "m13" "m17" "m19"
## [1] "m1"  "m5"  "m2"  "m10" "m7" 
## [1] "m1"  "m3"  "m20" "m15" "m8"

## [1] "K=3"
## [1] "m1"  "m5"  "m3"  "m16" "m11" "m13" "m17" "m9" 
## [1] "m1"  "m5"  "m3"  "m16" "m11" "m13" "m17" "m19"
## [1] "m1"  "m5"  "m2"  "m7"  "m10" "m4" 
## [1] "m1"  "m5"  "m3"  "m16" "m11" "m6"

## [1] "K=4"
## [1] "m1"  "m5"  "m3"  "m16" "m11" "m13" "m17" "m9" 
## [1] "m1"  "m5"  "m3"  "m16" "m11" "m13" "m17" "m19"
## [1] "m1"  "m5"  "m3"  "m16" "m11" "m6" 
## [1] "m1"  "m5"  "m3"  "m20" "m15" "m8"

Re-fitting zinbwave

Unsupervised

mod <- model.matrix(~ batch[our_cl != "-1"])
fn <- '../data/refit_zinbwave_slingshot.rda'
if (runZinb & !file.exists(fn)){
  zinbList <- lapply(Kvec, function(k){
    zinbFit(core[, our_cl != "-1"], X = mod, K = k, ncores = NCORES)
  })
  save(zinbList, file = fn)
}else{
  load(fn)
}

for(k in Kvec) {
  X <- getW(zinbList[[k - 1]])[, 1:k]

  lineages <- get_lineages(X, clus.labels = cl, start.clus = "m1")
  curves <- get_curves(X, clus.labels = cl, lineages = lineages)
  plot_curves(X, cl, curves, col.clus = pal)
  plot_tree(X, cl, lineages, col.clus = pal)

  print(paste0("K=", k))
  print(lineages$lineage1)
  print(lineages$lineage2)
  print(lineages$lineage3)
  print(lineages$lineage4)
}

## [1] "K=2"
##  [1] "m1"  "m2"  "m4"  "m16" "m12" "m14" "m19" "m17" "m9"  "m8" 
## [1] "m1"  "m2"  "m4"  "m16" "m6" 
## [1] "m1"  "m2"  "m4"  "m16" "m11"
## [1] "m1"  "m2"  "m4"  "m16" "m13"

## [1] "K=3"
##  [1] "m1"  "m2"  "m10" "m4"  "m16" "m14" "m19" "m17" "m9"  "m8" 
## [1] "m1"  "m2"  "m10" "m4"  "m16" "m11" "m6" 
## [1] "m1"  "m2"  "m10" "m4"  "m16" "m14" "m12"
## [1] "m1"  "m2"  "m10" "m4"  "m16" "m11" "m13"

## [1] "K=4"
##  [1] "m1"  "m2"  "m11" "m16" "m12" "m14" "m19" "m17" "m9"  "m8" 
## [1] "m1"  "m2"  "m11" "m16" "m4"  "m10"
## [1] "m1"  "m5"  "m3"  "m20" "m15"
## [1] "m1"  "m2"  "m11" "m6"

Supervised

for(k in Kvec){
  X <- getW(zinbList[[k - 1]])[, 1:k]

  lineages <- get_lineages(X, clus.labels = cl, start.clus = "m1",
                           end.clus = c("m8", "m4"))
  curves <- get_curves(X, clus.labels = cl, lineages = lineages)
  plot_curves(X, cl, curves, col.clus = pal)
  plot_tree(X, cl, lineages, col.clus = pal)

  print(paste0("K=", k))
  print(lineages$lineage1)
  print(lineages$lineage2)
  print(lineages$lineage3)
  print(lineages$lineage4)
}

## [1] "K=2"
## [1] "m1"  "m2"  "m16" "m12" "m14" "m19" "m17" "m9"  "m8" 
## [1] "m1"  "m5"  "m3"  "m15" "m20"
## [1] "m1"  "m2"  "m10" "m4" 
## [1] "m1"  "m2"  "m16" "m6"

## [1] "K=3"
## [1] "m1"  "m2"  "m11" "m16" "m14" "m19" "m17" "m9"  "m8" 
## [1] "m1"  "m2"  "m11" "m16" "m14" "m12"
## [1] "m1"  "m5"  "m3"  "m20" "m15"
## [1] "m1"  "m2"  "m11" "m16" "m18"

## [1] "K=4"
##  [1] "m1"  "m2"  "m11" "m16" "m12" "m14" "m19" "m17" "m9"  "m8" 
## [1] "m1"  "m2"  "m11" "m16" "m4" 
## [1] "m1"  "m5"  "m3"  "m20" "m15"
## [1] "m1"  "m2"  "m11" "m6"

DE analysis

Here is the kind of plots we want to be able to have

de <- read.csv('../data/oe_markers.txt', 
              stringsAsFactors = F, header = F)
de <- de$V1
nv <- assays(se)$normalizedValues
plotHeatmap(nv[rownames(nv) %in% de, ],
            main = 'Normalized values, 1000 most variable genes',
            clusterSamplesData = ceObj@dendro_samples,
            sampleData = as.matrix(sampleData$ours))

## Warning in .Method(..., deparse.level = deparse.level): number of rows of
## result is not a multiple of vector length (arg 1)

Session Info

sessionInfo()

## R version 3.3.3 (2017-03-06)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: CentOS release 6.9 (Final)
## 
## locale:
##  [1] LC_CTYPE=ja_JP.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=ja_JP.UTF-8        LC_COLLATE=ja_JP.UTF-8    
##  [5] LC_MONETARY=ja_JP.UTF-8    LC_MESSAGES=ja_JP.UTF-8   
##  [7] LC_PAPER=ja_JP.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=ja_JP.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    methods   stats     graphics  grDevices utils    
## [8] datasets  base     
## 
## other attached packages:
##  [1] rARPACK_0.11-0               digest_0.6.12               
##  [3] RColorBrewer_1.1-2           Rtsne_0.13                  
##  [5] matrixStats_0.52.2           magrittr_1.5                
##  [7] gplots_3.0.1                 ggplot2_2.2.1               
##  [9] slingshot_0.0.3-5            princurve_1.1-12            
## [11] zinbwave_0.99.1              clusterExperiment_1.3.0-9006
## [13] scone_0.99.0                 scRNAseq_1.0.0              
## [15] SummarizedExperiment_1.4.0   Biobase_2.34.0              
## [17] GenomicRanges_1.26.4         GenomeInfoDb_1.10.3         
## [19] IRanges_2.8.2                S4Vectors_0.12.2            
## [21] BiocGenerics_0.20.0          knitr_1.16                  
## 
## loaded via a namespace (and not attached):
##   [1] shinydashboard_0.6.0     R.utils_2.5.0           
##   [3] RSQLite_1.1-2            AnnotationDbi_1.36.2    
##   [5] htmlwidgets_0.8          grid_3.3.3              
##   [7] trimcluster_0.1-2        BiocParallel_1.8.2      
##   [9] RNeXML_2.0.7             DESeq_1.26.0            
##  [11] munsell_0.4.3            codetools_0.2-15        
##  [13] statmod_1.4.29           scran_1.2.2             
##  [15] DT_0.2                   miniUI_0.1.1            
##  [17] colorspace_1.3-2         uuid_0.1-2              
##  [19] pspline_1.0-17           robustbase_0.92-7       
##  [21] NMF_0.20.6               tximport_1.2.0          
##  [23] hwriter_1.3.2            rhdf5_2.18.0            
##  [25] rprojroot_1.2            EDASeq_2.8.0            
##  [27] diptest_0.75-7           R6_2.2.1                
##  [29] doParallel_1.0.10        ggbeeswarm_0.5.3        
##  [31] taxize_0.8.4             locfit_1.5-9.1          
##  [33] flexmix_2.3-14           bitops_1.0-6            
##  [35] reshape_0.8.6            assertthat_0.2.0        
##  [37] scales_0.4.1             nnet_7.3-12             
##  [39] beeswarm_0.2.3           gtable_0.2.0            
##  [41] phylobase_0.8.4          RUVSeq_1.8.0            
##  [43] bold_0.4.0               rlang_0.1.1             
##  [45] genefilter_1.56.0        splines_3.3.3           
##  [47] rtracklayer_1.34.2       lazyeval_0.2.0          
##  [49] hexbin_1.27.1            yaml_2.1.14             
##  [51] reshape2_1.4.2           abind_1.4-5             
##  [53] GenomicFeatures_1.26.4   backports_1.1.0         
##  [55] httpuv_1.3.3             tools_3.3.3             
##  [57] gridBase_0.4-7           dynamicTreeCut_1.63-1   
##  [59] stabledist_0.7-1         Rcpp_0.12.11            
##  [61] plyr_1.8.4               visNetwork_1.0.3        
##  [63] progress_1.1.2           zlibbioc_1.20.0         
##  [65] purrr_0.2.2.2            RCurl_1.95-4.8          
##  [67] prettyunits_1.0.2        viridis_0.4.0           
##  [69] zoo_1.8-0                cluster_2.0.6           
##  [71] RSpectra_0.12-0          data.table_1.10.4       
##  [73] mvtnorm_1.0-6            whisker_0.3-2           
##  [75] gsl_1.9-10.3             aroma.light_3.4.0       
##  [77] mime_0.5                 evaluate_0.10           
##  [79] xtable_1.8-2             XML_3.98-1.7            
##  [81] mclust_5.3               gridExtra_2.2.1         
##  [83] biomaRt_2.30.0           scater_1.2.0            
##  [85] tibble_1.3.1             KernSmooth_2.23-15      
##  [87] R.oo_1.21.0              htmltools_0.3.6         
##  [89] segmented_0.5-2.0        pcaPP_1.9-61            
##  [91] tidyr_0.6.3              geneplotter_1.52.0      
##  [93] howmany_0.3-1            DBI_0.6-1               
##  [95] MASS_7.3-47              fpc_2.1-10              
##  [97] MAST_1.0.5               boot_1.3-19             
##  [99] ShortRead_1.32.1         Matrix_1.2-10           
## [101] ade4_1.7-6               R.methodsS3_1.7.1       
## [103] gdata_2.17.0             igraph_1.0.1            
## [105] rncl_0.8.2               GenomicAlignments_1.10.1
## [107] registry_0.3             numDeriv_2016.8-1       
## [109] locfdr_1.1-8             plotly_4.6.0            
## [111] xml2_1.1.1               foreach_1.4.3           
## [113] annotate_1.52.1          vipor_0.4.5             
## [115] rngtools_1.2.4           pkgmaker_0.22           
## [117] XVector_0.14.1           stringr_1.2.0           
## [119] copula_0.999-16          ADGofTest_0.3           
## [121] softImpute_1.4           Biostrings_2.42.1       
## [123] rmarkdown_1.5            dendextend_1.5.2        
## [125] edgeR_3.16.5             kernlab_0.9-25          
## [127] shiny_1.0.3              Rsamtools_1.26.2        
## [129] gtools_3.5.0             modeltools_0.2-21       
## [131] rjson_0.2.15             nlme_3.1-131            
## [133] jsonlite_1.4             viridisLite_0.2.0       
## [135] limma_3.30.13            lattice_0.20-35         
## [137] httr_1.2.1               DEoptimR_1.0-8          
## [139] survival_2.41-3          prabclus_2.2-6          
## [141] iterators_1.0.8          glmnet_2.0-10           
## [143] class_7.3-14             stringi_1.1.5           
## [145] mixtools_1.1.0           latticeExtra_0.6-28     
## [147] caTools_1.17.1           memoise_1.1.0           
## [149] dplyr_0.5.0              ape_4.1